Task16b - Parallel Coordinates - Overiev of simulations
Setup and Read Data
source("ams_initialize_script.R")
source("SCIM_calculation.R")
source("ivsc_2cmt_RR_V1.R")
dirs$rscript_name = "Task16b_Parallel_Coordinates_Soluble_2019-11-12_AFIRvsSCIM_T0assumpt.R"
dirs$filename_prefix= str_extract(dirs$rscript_name,"^Task\\d\\d\\w?_")
data_in = read.csv("results/Task15_2019-11-12_10e3.csv",stringsAsFactors = FALSE)
data_in$id = 1:1e4Compute various quantities for comparing AFIR and SCIM, theory and simulation
Using adhoc theory calculation for SCIM
#put data into categories ----
data = data_in %>%
mutate(SCIM_thy= SCIM_adhoc_thy) %>%
mutate(Cavgss = dose_nmol/(CL*tau),
Kss_TL = Kd_TL + keTL/kon_TL,
Kss_DT = Kd_DT + keDT/kon_DT,
TL0 = T0*L0/Kss_TL,
koff_TL = Kd_TL*kon_TL,
koff_DT = Kd_DT*kon_DT,
ksynT = T0*keT + keTL*TL0,
ksynL = L0*(kon_TL*T0 + keL) - koff_TL*TL0,
Lss = ksynL/keL) %>%
mutate(AFIR_SCIM_sqerr = (AFIR_thy - SCIM_sim)^2) %>%
mutate(AFIRthy_category = case_when(AFIR_thy < 0.05 ~ "AFIRthy < 5%",
AFIR_thy > 0.30 ~ "AFIRthy > 30%",
AFIR_thy >= 0.05 & AFIR_thy <= 0.30 ~ "5% <= AFIRthy <= 30%"),
AFIRsim_category = case_when(AFIR_sim < 0.05 ~ "AFIRsim < 5%",
AFIR_sim > 0.30 ~ "AFIRsim > 30%",
AFIR_sim >= 0.05 & AFIR_sim <= 0.30 ~ "5% <= AFIRsim <= 30%"),
SCIMthy_category = case_when(SCIM_thy < 0.05 ~ "SCIMthy < 5%",
SCIM_thy > 0.30 ~ "SCIMthy > 30%",
SCIM_thy >= 0.05 & SCIM_thy <= 0.30 ~ "5% <= SCIMthy <= 30%"),
SCIMsim_category = case_when(SCIM_sim < 0.05 ~ "SCIMsim < 5%",
SCIM_sim > 0.30 ~ "SCIMsim > 30%",
SCIM_sim >= 0.05 & SCIM_sim <= 0.30 ~ "5% <= SCIMsim <= 30%"),
AFIRthy_AFIRsim_category = paste0(AFIRthy_category, ", ", AFIRsim_category),
AFIRthy_SCIMsim_category = paste0(AFIRthy_category, ", ", SCIMsim_category),
AFIRsim_SCIMsim_category = paste0(AFIRsim_category, ", ", SCIMsim_category),
SCIMthy_SCIMsim_category = paste0(SCIMthy_category, ", ", SCIMsim_category),
error_category = case_when(AFIR_SCIM_sqerr < 0.1 ~ "low_error",
TRUE ~ "high_error"))
data = data %>%
arrange(AFIR_thy) %>%
mutate(AFIRthy_category = factor(AFIRthy_category, levels = unique(AFIRthy_category))) %>%
arrange(AFIR_sim) %>%
mutate(AFIRsim_category = factor(AFIRsim_category, levels = unique(AFIRsim_category))) %>%
arrange(SCIM_sim) %>%
mutate(SCIMsim_category = factor(SCIMsim_category, levels = unique(SCIMsim_category)))
#check the assumptions of the data ----
data = data %>%
mutate(Ttotss = T0*Tfold,
koff_DT = Kd_DT*kon_DT,
assumption_plateau = AFIR_thy < 0.30,
assumption_drug_gg_T0 = Cavgss > 10*Ttotss,
assumption_drug_gg_KssDT = Cavgss > 10*Kss_DT,
assumption_koffDT_gt_keT = koff_DT > keT,
assumption_koffTL_fast = koff_TL > 1/30,
assumption_Cavgss_gg_LssKssDT_KssTL = Cavgss > 10*Kss_DT*Lss/Kss_TL,
assumption_T0simple = T0/(ksynT/keT) > 0.5 & T0/(ksynT/keT) < 2, #the simple formula works for T0
assumption_Lconst = Lss/L0 > 0.5 & L0/Lss < 2, #then SCIM = AFIR
assumption_Tss_gt_Lss = Tss_sim > Lss_sim,
assumption_all = assumption_plateau &
assumption_drug_gg_T0 &
assumption_drug_gg_KssDT &
assumption_koffDT_gt_keT &
assumption_koffTL_fast &
assumption_Cavgss_gg_LssKssDT_KssTL &
assumption_T0simple &
assumption_Tss_gt_Lss &
assumption_Lconst,
assumption_category = "",
assumption_category = ifelse(assumption_plateau,"","!plateau, "),
assumption_category = ifelse(assumption_drug_gg_T0,"","drug !>> T0"))Put data into error categories and summarize
threshold = 0.1
data_errss = data_in %>%
filter(abs(TLss_frac_change)>=threshold)
print(paste0(nrow(data_errss)," of ", nrow(data_in), " : Number of rows with TLss_frac_change > 0.1"))## [1] "208 of 10000 : Number of rows with TLss_frac_change > 0.1"data_err0 = data_in %>%
filter(abs(TL0_05tau_frac_change)>=threshold)
print(paste0(nrow(data_err0)," of ", nrow(data_in), " : Number of rows with TL0_05tau_frac_change > 0.1"))## [1] "0 of 10000 : Number of rows with TL0_05tau_frac_change > 0.1"# error historgram ----
data_quick_summ = data %>%
select(id,AFIR_thy, SCIM_sim, AFIR_SCIM_sqerr, TLss_frac_change, TL0_05tau_frac_change) %>%
gather(key,value,-c(id)) %>%
mutate(category = case_when((value < threshold) ~ "keep_low",
((value >= threshold) & (key %in% c("AFIR_SCIM_sqerr","SCIM_sim"))) ~ "keep_high",
((value >= threshold) & (key %in% c("AFIR_thy"))) ~ "keep_high_AFIR",
TRUE ~ "remove_high_error"))
g = ggplot(data_quick_summ, aes(value, fill = category))
g = g + geom_histogram()
g = g + facet_wrap(~key, scales = "free")
g = g + scale_fill_manual(values = c(keep_low = "grey80",
keep_high = "grey50",
remove_high_error = "red",
keep_high_AFIR = "blue"))
g = g + xgx_scale_x_log10()
g = g + ggtitle("")
print(g)
#keep only the simulations with no issues
data_keep = data %>%
filter(TLss_frac_change < threshold,
TL0_05tau_frac_change < threshold)
#put simulations into different categories
data_summary = data_keep %>%
group_by(AFIRthy_SCIMsim_category) %>%
count() %>%
arrange(desc(n))
kable(data_summary)| AFIRthy_SCIMsim_category | n |
|---|---|
| AFIRthy < 5%, SCIMsim < 5% | 4270 |
| AFIRthy < 5%, SCIMsim > 30% | 1667 |
| 5% <= AFIRthy <= 30%, SCIMsim > 30% | 1275 |
| 5% <= AFIRthy <= 30%, 5% <= SCIMsim <= 30% | 997 |
| AFIRthy > 30%, SCIMsim > 30% | 916 |
| AFIRthy < 5%, 5% <= SCIMsim <= 30% | 759 |
AFIRsim vs SCIMsim : 3x3 plot colors —-
param2uniform = function(x) {(log(x) - log(min(x)))/(log(max(x))-log(min(x)))}
data_plot = data_keep %>%
mutate_at(vars(AFIR:kon_TL,dose_mpk), funs(tf=param2uniform(.))) %>%
select(id,contains("AFIR"),contains("SCIM"), T0_tf:kon_TL_tf, dose_mpk_tf, contains("assumption")) %>%
gather(param,param_value,-c(id, contains("AFIR"), contains("SCIM"), contains("assumption"))) %>%
mutate(param = str_replace(param,"_tf",""))
#sort by average param value in one category to help with visualization ----
data_summ = data_plot %>%
filter(AFIRthy_SCIMsim_category == "AFIRthy < 5%, SCIMsim > 30%") %>%
group_by(param,AFIRthy_SCIMsim_category) %>%
summarise(x = mean(param_value)) %>%
arrange(x) %>%
ungroup()
kable(data_summ)| param | AFIRthy_SCIMsim_category | x |
|---|---|---|
| Kd_TL | AFIRthy < 5%, SCIMsim > 30% | 0.3746478 |
| keT | AFIRthy < 5%, SCIMsim > 30% | 0.4261922 |
| Kd_DT | AFIRthy < 5%, SCIMsim > 30% | 0.4417119 |
| Lfold | AFIRthy < 5%, SCIMsim > 30% | 0.4664135 |
| dose_mpk | AFIRthy < 5%, SCIMsim > 30% | 0.4688410 |
| kon_DT | AFIRthy < 5%, SCIMsim > 30% | 0.5014002 |
| keL | AFIRthy < 5%, SCIMsim > 30% | 0.5076666 |
| keDT | AFIRthy < 5%, SCIMsim > 30% | 0.5169649 |
| kon_TL | AFIRthy < 5%, SCIMsim > 30% | 0.5849314 |
| L0 | AFIRthy < 5%, SCIMsim > 30% | 0.6468130 |
| T0 | AFIRthy < 5%, SCIMsim > 30% | 0.6970433 |
data_plot = data_plot %>%
mutate(param = factor(param,
levels = data_summ$param))
g = ggplot(data_plot, aes(x=param,y=param_value, group = id, color = assumption_all, alpha = assumption_all))
g = g + geom_line()
g = g + facet_grid(SCIMsim_category~AFIRthy_category,switch = "y")
g = g + theme(axis.text.x = element_text(angle = 45, hjust = 1))
g = g + labs(x = "Parameter", y = "Parameter Value")
g = g + guides(colour = guide_legend(override.aes = list(alpha = 1)))
g = g + scale_color_manual(values = c(`TRUE` = "blue", `FALSE` = "red"))
g = g + scale_alpha_manual(values = c(`TRUE` = .1, `FALSE` = 0.01))
g = xgx_save(7,7,dirs,"Parallel_Coord_Soluble_3x3_AFIRthy_AFIRsim","")
g1 = g
print(g)
Focus on when assumptions are true and SCIM > 30% while AFIRthy < 5%
#explore data data where all assumptions are true and still
#AFIRsim > 30% and AFIRthy < 5% ---- on look, there is lots of L0!!!
#focus on this plot
data_new = data_plot %>%
filter(SCIMsim_category == "SCIMsim > 30%",
AFIRthy_category == "AFIRthy < 5%",
assumption_all == TRUE)
g = g1
g = g %+% data_new
g = g + geom_line(alpha = 0.05)
g = xgx_save(5,5,dirs,"Parallel_Coord_Soluble_AFIRthy_lt_5_SCIMsim_ge_30","")
g2= g
#print(g)Identify patients where all assumptions true and theory vs sim disagree.
source("SCIM_calculation.R")
source("ivsc_2cmt_RR_V1.R")
model = ivsc_2cmt_RR_KdT0L0()
library(RxODE)
id = unique(sort(data_new$id))
print("these IDs, even with all the restrictions, AFIR and SCIM still don't match")## [1] "these IDs, even with all the restrictions, AFIR and SCIM still don't match"print(id)## [1] 6 267 1289 2519 3186 4526 4621 5355 5484 5838 6009 6268 6322 7077
## [15] 7593 8379 8462 9268 9302 9624for (id_plot in id[1:5]) {
g = g2
g = g + geom_line(data = filter(data_new,id==id_plot),
size = 2,
color = "black")
g = g + ggtitle(paste("id =",id_plot))
print(g)
filepref = paste0("Parallel_Coord_Soluble_AFIRthy_lt_5_SCIMsim_ge_30_",id_plot)
g = xgx_save(5,5,dirs,filepref,"")
#simulate a patient where theory and simulation disagree
param = data %>%
filter(id==id_plot)
assumptions = param %>%
select(contains("assumption")) %>%
t()
kable(assumptions)
tmax = 365 #days
tau = param$tau #days
dose_nmol = param$dose_nmol
compartment = 2
infusion = TRUE
nam = names(param)
param_as_double = param %>%
as.numeric() %>%
setNames(nam)
param_as_double = param_as_double[model$pin]
param_print = param_as_double %>%
t() %>%
as.data.frame() %>%
mutate(CL = signif(keD/V1,2),
id = id_plot) %>%
select(id, CL,T0,L0,Kd_DT,Kd_TL,kon_DT,kon_TL,keT,keL,keDT,keTL)
ev = eventTable(amount.units="nmol", time.units="days")
sample.points = c(seq(0, tmax, 0.1), 10^(-3:0)) # sample time, increment by 0.1
sample.points = sort(sample.points)
sample.points = unique(sample.points)
ev$add.sampling(sample.points)
if (infusion == FALSE) {
ev$add.dosing(dose=dose_nmol, start.time = tau, nbr.doses=floor(tmax/tau), dosing.interval=tau, dosing.to=compartment)
} else {
ev$add.dosing(dose=dose_nmol, start.time = tau, nbr.doses=floor(tmax/tau)+1, dosing.interval=tau, dosing.to=compartment, dur = tau)
}
sim = lumped.parameters.simulation(model, param_as_double, dose_nmol, tmax, tau, compartment, infusion)
thy = lumped.parameters.theory ( param_as_double, dose_nmol, tau, infusion)
sim_rename = sim
nam = names(sim_rename) %>%
str_replace_all("_sim$","")
names(sim_rename) = nam
sim_rename$type = "sim"
thy_rename = thy
nam = names(thy_rename) %>%
str_replace_all("_thy$","")
names(thy_rename) = nam
thy_rename$type = "thy"
compare = bind_rows(sim_rename,thy_rename) %>%
select(type,Dss,T0,L0,TL0,Ttotss,Lss,TLss,AFIR,SCIM)
init = model$init(param_as_double)
out = model$rxode$solve(model$repar(param_as_double), ev, init)
out = model$rxout(out)
out_plot = out %>%
select(time,D,T,DT,L,TL) %>%
gather(cmt,value,-time)
out_last = out_plot[(out$time==max(out$time)),]
g = ggplot(out_plot,aes(x=time,y=value, color = cmt, group= cmt))
g = g + geom_line()
g = g + geom_label(data = out_last, aes(label = cmt), show.legend = FALSE, hjust=1)
g = g + geom_vline(xintercept = tau, linetype = "dotted")
g = g + xgx_scale_x_time_units(units_dataset = "days", units_plot = "weeks")
g = g + xgx_scale_y_log10()
g = g + labs(y = "Concentration (nm)", color = "")
g = g + ggtitle(paste0( "id = ",param$id,
"\nAFIR_thy = ",signif(thy$AFIR_thy,2),
"\nAFIR_sim = ",signif(sim$AFIR_sim,2),
"\nSCIM_thy = ",signif(thy$SCIM_adhoc_thy,2),
"\nSCIM_sim = ",signif(sim$SCIM_sim,2)))
filepref = paste0("Parallel_Coord_Soluble_AFIRthy_lt_5_SCIMsim_ge_30_",id_plot)
g = xgx_save(5,5,dirs,filepref,"")
print(g)
#unfortunately, kable does not work properly inside for loop
print(t(param_print))
print(t(compare))
}
## [,1]
## id 6.00e+00
## CL 2.70e-01
## T0 2.58e+01
## L0 1.57e-02
## Kd_DT 1.61e+01
## Kd_TL 1.04e-01
## kon_DT 3.83e+02
## kon_TL 2.82e+01
## keT 1.32e+03
## keL 4.40e+00
## keDT 7.06e-02
## keTL 1.27e+00
## [,1] [,2]
## type "sim" "thy"
## Dss "6645781" "6686508"
## T0 "25.8" "25.8"
## L0 "0.0157" "0.0157"
## TL0 "2.717877" "2.717877"
## Ttotss "461494.0" "482428.5"
## Lss "0.2528003" "0.8001780"
## TLss " 1.896427" "-99.000000"
## AFIR "0.04333395" "0.04308006"
## SCIM "0.6977604" "2.2946954"
## [,1]
## id 2.67e+02
## CL 2.70e-01
## T0 1.61e+01
## L0 5.44e-03
## Kd_DT 8.58e+00
## Kd_TL 6.55e-01
## kon_DT 1.53e+03
## kon_TL 5.60e+00
## keT 1.17e+03
## keL 5.92e+00
## keDT 8.16e-02
## keTL 4.76e+00
## [,1] [,2]
## type "sim" "thy"
## Dss "2457735" "2480159"
## T0 "16.1" "16.1"
## L0 "0.00544" "0.00544"
## TL0 "0.05819535" "0.05819535"
## Ttotss "219844.3" "230849.0"
## Lss "0.03704329" "0.05223221"
## TLss " 0.01889042" "-99.00000000"
## AFIR "0.04766974" "0.04725864"
## SCIM "0.3246037" "0.4762670"
## [,1]
## id 1.289e+03
## CL 2.700e-01
## T0 1.410e+02
## L0 9.100e-04
## Kd_DT 6.030e+01
## Kd_TL 1.040e-02
## kon_DT 2.160e+01
## kon_TL 3.500e+03
## keT 1.530e+00
## keL 3.130e+00
## keDT 4.520e-02
## keTL 3.770e-01
## [,1] [,2]
## type "sim" "thy"
## Dss "319176.0" "319444.4"
## T0 "141" "141"
## L0 "0.00091" "0.00091"
## TL0 "12.21103" "12.21103"
## Ttotss "4763.398" "4874.636"
## Lss "0.1303555" "1.4716952"
## TLss " 11.13632" "-99.00000"
## AFIR "0.006366509" "0.006349270"
## SCIM " 0.9119888" "10.2807213"
## [,1]
## id 2.519e+03
## CL 2.700e-01
## T0 1.990e+01
## L0 9.150e-04
## Kd_DT 1.420e+01
## Kd_TL 2.160e-03
## kon_DT 5.690e+02
## kon_TL 6.360e+02
## keT 2.350e+02
## keL 4.170e+01
## keDT 1.020e-01
## keTL 4.480e+01
## [,1] [,2]
## type "sim" "thy"
## Dss "5867197" "5873016"
## T0 "19.9" "19.9"
## L0 "0.000915" "0.000915"
## TL0 "0.2508049" "0.2508049"
## Ttotss "45635.42" "45958.20"
## Lss "0.1026274" "0.2703649"
## TLss " 0.1561306" "-99.0000000"
## AFIR "0.005550212" "0.005539720"
## SCIM "0.6225181" "1.6499348"
## [,1]
## id 3.186e+03
## CL 2.700e-01
## T0 5.680e+01
## L0 1.550e-03
## Kd_DT 7.720e+01
## Kd_TL 1.050e-01
## kon_DT 5.080e+00
## kon_TL 1.760e+01
## keT 2.540e+01
## keL 1.980e+00
## keDT 6.660e-02
## keTL 2.660e-01
## [,1] [,2]
## type "sim" "thy"
## Dss "821671.6" "823412.7"
## T0 "56.8" "56.8"
## L0 "0.00155" "0.00155"
## TL0 "0.7329726" "0.7329726"
## Ttotss "20916.30" "21665.39"
## Lss "0.03127469" "0.10002005"
## TLss " 0.5117136" "-99.0000000"
## AFIR "0.03460014" "0.03452809"
## SCIM "0.6981347" "2.3076631"